
;;================================================================
;; Routine to generate the ATLASGAL - BGPS latitude comparison for the
;;   overlap between the two surveys.  BOTH catalogs have been
;;   restricted to the region -10 < l < 21 (something not done in
;;   Fig. 9 of {Contreras:2013}).


;; Open file for plotting
myps,'./emaf_paper/plots/atlasgal_glat.eps',/bw

deg = cgSymbol('deg')

;; Why use a stolen version of the plot, when you have the real data???
;; readcol,'./ancillary/ATLASGAL_GLAT.txt',b,n
csc = read_mrt('./ancillary/ATLASGAL_csc.dat')

;; Cut to BGPS overlap
ind = where(csc.glon GE 100)
csc[ind].glon -= 360.
ind = where(csc.glon GE -10.5)
csc = csc[ind]

;; Plot!
plothist,csc.glat,b,n,bin=0.1,thick=5,xtit='GLAT [deg]',/fill,fcolor='BLK4',$
         ytit='N per 0.1'+deg+' bin',charsize=1.,yr=[-100,600],/yst,color='BLK4'

;; Horizontal dashed line at N=0
;; vline,0,thick=3,/h


;; Make Gaussian fit for ATLASGAL data
yf = mpfitpeak(b,n,A,nt=3)
xc = findgen(4001)/1000. - 2.
yc = gauss_1(xc,A)
cgOplot,xc,yc,thick=5,color='BLK6',linestyle=2

;; ATLASGAL Legend
al_legend,/top,/left,box=0,textcolor=['Opposite','BLK4','BLK6'],$
          ['ATLASGAL','CENTROID = '+string(A[1],format="(F0.2)")+deg,$
           'FWHM = '+string(A[2]*2.355,format="(F0.2)")+deg],linsize=0.5,$
          color=['Background','BLK4','BLK6'],thick=5,linestyle=[0,0,2]


;;====================================
;; Moving on to the BGPS...
s = read_bgps_csv()
ind = where(s.glon_peak GE 250)
s[ind].glon_peak -= 360.
ind = where(s.glon_peak LE 21)
s = s[ind]

;; Make histogram data, then re-scale to match ATLASGAL
plothist,s.glat_peak,sb,sn,/noplot,bin=0.1
;; sn = sn*max(n)/max(sn)
cgOplot,sb,sn,thick=5,color='Opposite',psym=10

;; BGPS Gaussian
yf = mpfitpeak(sb,sn,A,nt=3)
yc = gauss_1(xc,A)
cgOplot,xc,yc,thick=5,color='Opposite'

;; BGPS Legend
al_legend,/top,/right,box=0,$
          textcolor=['Opposite','Opposite','Opposite'],$
          ['BGPS','CENTROID = '+string(A[1],format="(F0.2)")+deg,$
           'FWHM = '+string(A[2]*2.355,format="(F0.2)")+deg],linsize=0.5,$
          color=['Background','Opposite','Opposite'],thick=5,linestyle=[0,0,0]
          

;; Now, work the magic to create the difference plot...
sin = interpol(sn,sb,b)
sin[where(b LT -1. OR b GT 1.1)] = 0

diff = n - sin
print,m4_stat(diff),total(diff)/total(n)

;; Plot!
cgoplot,b,diff,color='Opposite',thick=5,psym=10,linestyle=5


al_legend,/bottom,/right,linesty=5,color='Opposite',thick=5,box=0,$
          ['ATLASGAL-BGPS'],charsize=0.8

al_legend,/bottom,/left,box=0,charsize=0.8,$
          ['(A-B)/A = '+string(total(diff)/total(n),format="(F0.2)")]

myps,/done

;;===============================================================
;; L-B plot

myps,'./emaf_paper/plots/atlasgal_scatter.eps',xsize=9,ysize=6.5,ct=41

bcolor='blue violet'

range=[-0.1,2.4]

wco = alog10(readfits('./local/co_surveys/Wco_DHT2001.fits',cohdr))

plotmap,wco,cohdr,xrange=[22,-11.5],yrange=[-1.75,1.75],charsize=1.1,$
        range=range,xmargin=[10,10],axsel=3,min_dpi=300,/notsquare


cgOplot,csc.glon,csc.glat,psym=16,symsize=0.4

cgOplot,s.glon_peak,s.glat_peak,psym=16,symsize=0.4,color=bcolor

integral = '!9' + string("362B) + '!X'
;;"
cgLoadct,41,/silent
cgcolorbar,/vertical,charsize=1.0,/right,range=range,$
           position=[0.915,!y.window[0]+0.1,0.945,!y.window[1]-0.1],minor=5,$
           title='log '+integral+' T!dmb!ndv [K km s!u-1!n]',$
           xthick=3,ythick=3,divis=0

al_legend,/top,/left,background_color='Background',psym=16,symsize=2.0*0.4,$
          color=['Black',bcolor],['ATLASGAL','BGPS']

;; bpa = fltarr(17)
;; bpb = fltarr(17)
;; i=0
;; ell = dindgen(17)*2.-11.
;; print,m4_stat(ell)

;; FOR i=0,n_elements(bp)-1 DO BEGIN
   
;;    ind = where( abs(csc.glon-ell[i]) LE 1.0, nind)
;;    bpa[i] = median(csc[ind].glat)
   
;;    ind = where( abs(s.glon_peak-ell[i]) LE 1.0, nind)
;;    bpb[i] = median(s[ind].glat_peak)
   
;;    i++
;; ENDFOR

;; cgOplot,ell,bpa,psym=10,thick=5
;; cgOplot,ell,bpb,psym=10,thick=5,color='dodger blue'



myps,/done

END
